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(57) Abstract: A system and a method solve the estimation problem of finding reflectance and illumination L. The system and 
method solve a functional of the unknown illumination L such that a minimum of the functional is assumed to yield a good estimate 
— ^ of the illumination Having a good estimate of the illumination L implies a good estimate of the reflectance R. The functional uses 
2 a variational framework to express requirements for the optimal solution. The requirements include: 1) that the illumination L is 
spatially smooth; 2) that the reflectance values are in the interval [0,1]- thus, when decomposing the image S, the solution should 
^ satisfy the constraint L>;3) that among all possible solutions, the estimate of the illumination L should be as close as possible to the 
^ image 5, so that the contrast of the obtained R is maximal; and 4) that the reflectance R complies with typical natural image behavior 
^ (e.g., the reflectance is piece -wise smooth). 



WO 02/089062 



PCT/US02/14987 



1 System and Method for Image Enhancement, Dynamic Range Compensation, 

2 and Illumination Correction 

3 Technical Field 

4 The technical field is enhancement of digital images. 

5 Background 

6 The human visual system is typically able to adapt to lighting variations across 



7 scenes, perceiving details in regions with very different dynamic ranges. Most image 

8 recording systems, however, fail this dynamic range compression task. As a result, images 

9 produced by these image recording systems are often of poor quality, compared to images 

10 produced by human perception. Another task that is often poorly performed by the image 

1 1 recording systems is that of color constancy. Humans perceive color in a way that is fairly 

12 independent of scene illumination, whereas the image recording systems are strongly 

13 influenced by spectral shifts. 



14 The above problems can be stated mathematically by describing a relationship 

15 between an acquired image S, a reflectance of objects of the image R, and an illumination L 

16 in a pixel-wise multiplication, or: 

17 S = LR. 

18 This expression means that at each point in the imageS, the color value is the multiplication 

19 of the reflectance value by the illumination value. Given an image S, the problem to be solved 

20 is removal of the effects of illumination and recovery of the reflectance image. That is, given 

21 S, find both R and L. However, calculation of both R and L is typically not possible. The 

22 solution then, is to generate a methodology that can estimate R and L 

23 Retinex theory deals with compensation for illumination in images. The goal of the 



24 Retinex theory is to decompose a given image S into the reflectance image/?, and the 

25 illumination image Z, such that at each point (x,y) in the image domain, S(x,y) equals the 

26 product of R(x,y) and L(x,y). The benefits of such decomposition include the possibility of 

27 removing illumination effects of back/front lighting, enhancing photographs and other image 

28 capture methods that include spatially varying illumination, such as images that contain 

29 indoor and outdoor zones, and correcting colors in images by removing illumination-induced 

30 color shifts. 

3 1 As noted above, recovering the illumination L from a given image S is known to be a 

32 mathematically ill-posed problem, and known algorithms vary in the manner and 

33 effectiveness of overcoming this limitation. The Retinex approach provides the framework 

34 for one such method. The Retinex methodology was motivated by Edward Land's research 
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1 of the human visual system, which is described in R.H. Land, "The Retinex Theory of Color 

2 Vision," ScLAmer., Vol. 237, pp. 108-128(1977). 

3 The first Retinex algorithms were of the random walk type. Subsequent Retinex 

4 algorithms used homomorphic filtering. Yet another group of Retinex algorithms are based 

5 on solving a Poisson equation. 

6 Figure 1 is a block diagram of an algorithm 10 that is representative of the general 

7 class of prior-art Retinex algorithms. In Figure 1, an input image S 1 1 is applied to a 

8 logarithmetic module 12 to produce a logarithmetic version 13 of the input imaged and its 

9 two components, illumination L and reflectance R. That is, s = log S> I = log L 9 and r = log R 9 

10 and thereby, s - / + r. Using the logarithmetic version 13 of the input images, an estimator 

1 1 module 14 computes an estimate 17 of the illumination, designated in Figure 1 as/*. The 

12 estimate 17 (/*) is then combined with the logarithmetic version 13 (s) at summer 18 to 

13 produce an estimate 19 of the reflectance (designated r*). Finally, the estimate 19 fr*) is 

14 converted from a logarithm to a number (antilogarithm) corresponding to the logarithm at 

15 expander 20, to produce a real number value as an estimate 21 of the reflectance (designated 

16 as R*). Prior art Retinex algorithms usually employ the same process as shown in Figure 1 . 

1 7 The Retinex-based algorithms take several different forms. One such form is the 

1 8 random walk algorithm, which is a discrete time random process in which the "next pixel 

19 position" is chosen randomly from neighbors of the current pixel position. Random walk 

20 type Retinex algorithms are variants of the following basic formulation: A large number of 

21 walkers are initiated at random locations of the logarithmetic version 13 (s), adopting a gray- 

22 value of their initial position. An accumulator imaged that has the same size as s is 

23 initialized to zero. As a walker walks around, the walker updates A by adding the value of 

24 the walker to each position (x t y) that the walker visits. The illumination image is obtained by 

25 normalizing the accumulator imaged, i.e., the value at each position of the accumulator 

26 image A is divided by the number of walkers that visited that position. 

27 By using many walkers with long paths, one can easily verify that the accumulator 

28 value at each position converges to a Gaussian average of that position's neighbors. 

29 Another type of Retinex algorithm uses homomorphic filtering, where a low-pass 

30 filter is used to reconstruct / from s. Homomorphic Retinex algorithms are based on the fact 

3 1 that the reflectance image R corresponds to the sharp details in the image (i.e., the edges), 

32 whereas the illumination L is expected to be spatially smooth. Then, a reasonable guess for / 

33 is /* = LP{s} 9 where LP is a convolution with a wide Gaussian kernel. Thus, the Retinex 
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1 algorithm using homomorphic filtering actually applies the same process as the random walk 

2 algorithms by a single direct convolution. 

3 Since the illumination L is expected to be spatially smooth, the derivative of the 

4 illumination L should be close to zero. However, considering the assumption that the 

5 reflectance R is piece-wise constant, the derivative of the reflectance R should vanish almost 

6 everywhere, with high values along edges of an image. Taking the derivative of the sum s= I 

7 + r and clipping out high derivative peaks, implies that the clipped derivative signal 

8 corresponds only to the illumination L. Poisson equation-type Retinex algorithms that rely on 

9 the Mondrian world model, use the above assumptions on the reflectance Ras a piece-wise 
10 constant image. Applying the Laplacian, and the following clipping operation: 

11 

12 t(As) = As where \As\ <t 

13 x(As) = 0 otherwise, 

14 yields the following Poisson equation 

15 A/* = t(Ay) 

16 A solution to the Poisson equation may involve an iterative procedure that effectively inverts 

17 the Laplacian operator. Improvements to the method involve extracting discontinuities from 

1 8 the image gradient magnitude instead of the Laplacian to provide better boundary conditions. 

19 One solution involves use of an iterative algorithm having a "reset" non-linearity that 

20 enforces the constraint / >s. The algorithm performs the interactive procedure. 

2 * 2 J 

22 where D n is a translation operator, shifting the image s by the rf 1 element of a sequence of 

23 spirally decaying translation vectors. Removing the max operation yields a simplified 

24 version 

" 2 

26 The above equation is a simple averaging operation that smoothes the images. The non- 
27 linear (max) operation forces the illumination image / to satisfy the constraint/* 2> s. 

28 Despite use of these algorithms, current image recoding and image enhancement 

29 systems and methods cannot produce images that are of sufficient quality to be comparable to 

30 images as perceived by the human vision system. Thus, an improved method and ystem are 

3 1 required to better remove illumination effects and to recover the reflectance image. 
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1 Summary 

2 A system and a method solve the estimation problem of finding reflectance/? and 

3 illumination L. The system and method use a functional of the unknown illumination L such 

4 that a minimum of the functional is assumed to yield a good estimate of the illumination!. 

5 Having a good estimate of the illumination! implies a good estimate of the reflectance R. 

6 The functional uses a variational framework to express requirements for the optimal 



7 solution. The requirements include: 1) that the illumination L is spatially smooth; 2) that the 

8 reflectance values are in the interval [0,1] - thus, when decomposing the image S, the solution 

9 should satisfy the constraint L> S; 3) that among all possible solutions, the estimate of the 

1 0 illumination L should be as close as possible to the imaged, so that the contrast of the 

1 1 obtained R is maximal; and 4) that the reflectance R complies with typical natural image 

12 behavior (e.g., the reflectance is piece-wise smooth). 



13 The four requirements lead to a well-behaved optimization problem having quadratic 

14 programming, which is convex, with a single solution. The method and the system use a 

1 5 numerical algorithm for the solution, where the algorithm is based on a multi-resolution 

1 6 decomposition of the image S. 

17 The method and the system are adaptable to different image situations. In a first 



18 situation, given a color image in the RGB color-space, applying the algorithm on the three 

19 color channels separately produces color correction. An alternative embodiment involves 

20 converting the input image S into a chromatic-luminance color space. The algorithm is then 

21 applied to the luminance alone, resulting in an image that preserves the input colors, and 

22 improves the local contrast. 

23 A second situation accounts for the fact that although the illumination L is assumed to 

24 be spatially smooth on a global level, the assumption is sometimes wrong. This situation is 

25 shown dramatically in Figure 9, where the illumination L is expected to have a sharp edge on 

26 the window boundaries. In situations such as that shown by the example of Figure 9, 

27 allowing the illumination L to include sharp edges may be enabled by forcing piece-wise 

28 smoothness on the illumination L, rather than assuming global smoothness. 



29 A third situation does not use the assumption of piece-wise smoothness for the 

30 reflectance R. Instead, a better quality prior imageS is used. The better prior S allows use of 

3 1 multi-resolution and unisotropic considerations. 

32 A fourth and final situation accounts for the fact that some illumination effect is 

33 needed in order to give the human user a natural feel of the recorded image. Thus, the 
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1 method and the system may include a process that partiallyreturns the illumination such that 

2 a final output image includes some illumination, but with improved contrast. 

3 Using the four requirements enumerated above, the variational framework and the 

4 algorithm provide a novel solution to the problem of estimating illumination. The variational 

5 framework and the algorithm may be applied in many image recording situations, including 

6 in digital cameras as an automatic illumination compensation module, aimed at improving 

7 image reproduction quality, and in scanners and printers, as a special effect that improves the 

8 visual quality of a scanned/printed image. 

9 Description of the Drawings 

10 The detailed description will refer to the following drawings, in which like numbers 

1 1 refer to like objects, and in which: 

1 2 Figure 1 is a diagram of a prior art algorithm for estimating illumination; 

1 3 Figure 2 is a block diagram of an improved algorithm for decomposing an image; 

14 Figure 3 is a diagram of a system that uses the improved algorithm of Figure 2; 

1 5 Figure 4 is a flowchart showing a process according to the algorithm of Figure 2; 

1 6 Figure 5 is a flowchart of a sub-routine of the process of Figure 4; 

1 7 Figure 6 is a flowchart of a further sub-routine of Figure 4; 

1 8 Figure 7 is a block diagram of an alternate improved algorithm; 

19 Figure 8 is a block diagram of yet another alternate improved algorithm; and 

20 Figure 9 is an image triplet showing the effects of illumination on an image. 

21 Detailed Description 

22 The human visual system is typically able to adapt to lighting variations acrossscenes 



23 (for example, shadows, strong illumination source), perceiving details in regions with very 

24 different dynamic ranges. Most image recording systems (e.g., cameras), however, fail this 

25 dynamic range compression task. As a result, images produced by these image recording 

26 systems are often of poor quality, compared to images produced by human perception. 

27 Another task that is often poorly performed by the image recording systems is that of color 

28 constancy. Humans perceive color in a way that is fairly independent of scene illumination, 

29 whereas the images recording systems are strongly influenced by spectral shifts. 



30 The above problems can be mathematically formulated by describing a relationship 

31 between an acquired image S 9 a reflectance R of objectsof the image, and an illumination L 

32 in a pixel-wise multiplication, or: 

33 S = LR. 
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1 This expression means that at each point in the imaged, the color value is the multiplication 

2 of the reflectance value by the illumination value. The effect of this relationship can be seen 

3 by a simple experiment. A person or other object in a room is placed near a window (or other 

4 light source) during daylight such that a strong light source (i.e., bright daylight) illuminates 

5 part of an image S to be recorded. Thus, the illumination is strong through the window and 

6 weak inside the room. A photograph is then taken, and the recorded imageS is decomposed 

7 into reflectance R and illumination L. The person can barely be seen in the image S and the 

8 illumination L, but can be seen in the reflectance R. 

9 Given an image 5, the problem to be solved is removal of the effects of illumination 

10 and recovery of the reflectance image. That is, givenS, find both R and L. However, 

1 1 calculation of either R or L is typically not possible. The solution then, is to generate a 

12 methodology that can estimate R and L. 

13 Figure 2 is a block diagram of an improved algorithm 1 00 that may be used to 

14 estimate the illumination L and the reflectance R. In Figure 2, an image S 101 is input to a 

15 logarithm module 102 to produce a logarithm s 103 of the image 5101 and its two 

1 6 components, illumination I and reflectance R. That is, s = log S, I = log L 9 and r = log R 9 and 

17 thereby, s = / + r. Using the logarithm s 103, an image processing module 104 uses a 

18 Projected Normal Steepest Descent or similar algorithm, with multi-resolution processing, to 

19 compute an estimate 107 of the illumination, designated in Figure 2 as/*. The estimate 107 

20 (/*) is then combined with the output 103 (?) at summer 108 to produce an estimate 109 of 

21 the reflectance (designated r*). Finally, the estimate 109 (r*) is converted from a logarithm 

22 to its corresponding base number value at exponential module 1 10, to produce a number 

23 value as an estimate 1 1 1 of the reflectance (designated as R*). 

24 Figure 3 is a block diagram of a representative component (a digital camera) 1 1 2 that 

25 uses the system and the method described herein. The component 1 12 includes an image 

26 capture device 1 14 that captures an image. The image capture device may include lenses, a 

27 CCD array, and other sub-components. An image converter 115 may be an optional 

28 subcomponent used for pre-processing of the captured image. An image processor 116 may 

29 include the necessary processing software and algorithms to enhance the captured image 

30 according to the methods and systems described herein. The component 1 12 may include an 

31 image memory 1 17 that stores captured images, both before and after processing by the 

32 image processor 1 16. Finally, a display device 118, which may be a liquid crystal display 

33 device, for example, provides a visual display for a user of the component 1 12. 
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1 To arrive at the solution to the problem of estimating/? and L, a method and a system 

2 for applying the method begins by formulating the estimation problem of finding R and L 

3 from S as an optimization problem. That is, a functional of one of the unknowns (either lor 

4 R) is provided such that a minimum of the functional yields the desired result. The 

5 formulation uses a variational framework to express requirements from the optimal solution. 

6 The framework embodies the following requirements: 

7 1. The illumination L is spatially smooth. 

8 2. The reflectance R is restricted to the unit interval, which adds the constraint L> S. 

9 Since the log function is monotone, / > s. 

10 3. Setting / = Const, where Const is any constant above the maximal value ofr, yields a 

1 1 trivial solution that satisfies the two previous requirements. Thus, a requirement is 

12 added such that the illumination image / is close to the intensity image s y i.e., the 

13 value of / minimizes a penalty term of the form dist s), e.g., the£2 norm (/ - s)2. 

14 4. The reflectance image s = / - r can be assumed to have a high prior probability. One 

15 of the simplest prior functions used for natural images assigns high probability to 

16 spatially smooth images. 

17 5. The illumination continues smoothly as a constant beyond the image boundaries. 

18 This is an artificial assumption required for boundary conditions that would have 

1 9 minor effect on the final results. 

20 Collecting all the above requirements into one expression yields the following penalty 

21 functional: 

Minimize: p ^ = jj v/ |2 +a (j_ s y +>? | V (/- ^)| 2 )&^ (1) 
Subjectto: ,^ \ nd ^^ 0 on 3ft, 

22 where ft is the support of the image, 3 ft is the image boundary, and n is the normal to the 

23 boundary, a and P are free non-negative real parameters. In the functional F[/), the first 

24 penalty term |v/| 2 ) forces spatial smoothness on the illumination image. This choice of 

25 smoothness penalty is natural, keeping in mind that minimizing j|v/| 2 }&rfy translates into 

26 the Euler-Lagrange (EL) equation A/ = 0. The steepest descent solution to the first penalty 

27 term is a Gaussian smoothing operation with increasing variance of the initial condition. As 

28 mentioned in the previous section, several authors proposed Gaussian smoothing of jfor the 

29 illumination reconstruction. 
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1 The second penalty term (/ - s) 2 forces a proximity between / and s. The difference 

2 between these images is exactly r, which means that the norm ofr should be small (i.e., R 

3 tends to white). This term is weighted by the free parameter a The main objective of this 

4 term is a regularization of the problem that makes the second penalty term better conditioned. 

5 Notice that, in addition, the solution / is forced to be / > s. 

6 The third penalty term is motivated by the Bayesian prior probability theory and 

7 encourages the reflectance image r to be a "visually pleasing image." In practice, this prior 

8 probability expression forces r to be spatially smooth. The third penalty term is weighted by 

9 the free parameter p. More complicated Bayesian expressions may be used allowing sharp 

10 edges, textures, 1 / /behavior, etc. As long as this expression is purely quadratic, the above 

1 1 minimization problem remains fairly simple. 

12 The above-defined problem has a Quadratic Programming (QP) form. The necessary 

13 and sufficient conditions for the minimization problem are obtained with the Euler-Lagrange 

14 equations: 



dF ^ = 0 = - A/ + a(l -s)- J3A(l -s)andl>s 



dl 

or 

l=s 



(2) 



1 5 The differential equation does not have to hold when / = s. 

16 The minimization problem is QP with respect to the unknown image /. A Projected 

17 Normalized Steepest Descent (PNSD) algorithm, accelerated by a multi-resolution technique, 

18 is used to solve this minimization problem. 

19 The PNSD algorithm requires the application of a Normalized Steepest Descent 

20 (NSD) iteration that minimizes the functional F[/] 5 followed by a projection onto the 

21 constraints. A NSD iteration has the format: 

22 h m h-i~ MnsdG 

23 where lj and lj _ ] are the illumination images at step/* and j - 7, respectively, G is the 

24 gradient of F[l\ 9 and /i NSD ' s the optimal line-search step size. For Equation (2), the 

25 gradient of F[[\ is given by: 

26 G= V/;.; + ^a- fiA)(lj - . - 5), 
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K 

2 Using integration by parts, I I V G 1 2 = - J GAG up to the boundary conditions. 

3 An alternative approach is the Steepest Descent (SD) algorithm, where /j nsD is 

4 replaced by a constant value // $D> suc ^ t ^at: 



5 where ^{,4} refers to the greatest eigenvalue of the linear operator^. This alternative 

6 method saves computations at the expense of a slightly slower convergence. 

7 Finally, projecting onto the constraint / > s is done by /j = max(/,-, s). In practice, G 

8 can be calculated by: 

9 G = G A + OL(lj-l-s)-&(G A -G B ) 

10 where 

11 G A AA/y_;, 

12 Gfl AAjJt 

13 Similarly, /i is given by: 



14 where 

^ 4 J|G|\ 

15 The Laplacian may be approximated by a linear convolution with the kernel ki^p 
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1 and the integrations are approximated by summations: 



n m 

J|VG|* = - JCAG 



2 where G[m, n] = G(mAx, «Ay). In order to accommodate the boundary conditions, as given 

3 in Equation (1), the above convolution is applied on an expanded version of the image G 

4 The extension is done by replicating the first and last columns and rows. After convolution, 

5 the additional rows and columns are removed. 

6 Although simple, the PNSD algorithm usually converges slowly. Instead of general 

7 acceleration schemes, the fact that the unknown image / is assumed to be smooth may be 

8 used. Specifically, the method applies a multi-resolution algorithm that starts by estimating a 

9 coarse resolution image /, expands the resolution image / by interpolation, and uses the result 

10 as an initialization for the next resolution layer. This way, only a few iterations at each 

1 1 resolution are required for convergence. 

12 Summarizing the above, a block diagram of an algorithm 120 for the solution of 

13 Equation (1) is shown in Figure 4. The algorithm 120 begins with input block 121, where an 

14 input to the algorithm 120 is an images of size [N, M\ 9 and the two parameters a and p. 

15 In initialization block 123, a Gaussian pyramid of the image s is computed. The thus- 

16 constructed pyramid contains p resolution layers (from fine (1 ) to coarse (p)) with the current 

17 resolution layer (k) set to p. In block 125, Tk iterations of tte PNSD algorithm are applied to 

18 the kth resolution layer, until all resolution layers are checked, block 127. In block 129, the 

19 next resolution layer is updated. When all resolution layers are processed, the result is the 

20 final output of the algorithm 120. 

21 Figure 5 is a block diagram showing the initialization routine 123 in more detail. In 

22 block 131, a counter is initialized with 1=1. In block 132, the Gaussian pyramid is 

23 constructed by smoothing the image with a kernel, such as the kernel kpYR* 
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1 In block 1 33, the image is decimated by, for example, a 2: 1 ratio. The process is repeated 

2 (block 134) and the counter is incremented (block 135) p times where for this section p < Ig 

3 (min (M, N)). This process produces a sequence of images {s k conventionally named the 

4 "Gaussian Pyramid of s." The image sj is the original image s, and Sp is an image with the 

5 coarsest resolution for the Gaussian pyramid. 

6 In block 137, the process sets* =p, i.e., starts at the coarsest resolution layer k, and 

7 sets the initial condition Iq = max {sp} . 

8 Figure 6 is a block diagram of the PNSD processing main routine 125. The routine 

9 125 starts at block 141, where for the kfh resolution layer, the routine 125 calculates: 



G B =A k s k 

1 0 Where: A^= c, * kup2~^ & - 1\ namely a convolution of <; with the kernel hup as 

1 1 specified above, normalized by a factor 2 (2<l " ,)) . 



12 



13 



14 



Then, for/ = 1, ..,r k : 

In block 143, the routine 125 calculates the gradient: 

Ga = V/-.. 

G = G A +a (/,_, -s k )- fi(G A -G B ). 
In block 145, the routine 125 calculates n nsd : 



M* = (G,G), 

fi B - -{g,a k g), 



N M 



15 Where: 



16 



In block 147, the routine 125 completes the NSD iteration: 
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In block 149, the result is projected onto the constraints: 

/ y =max^t 

The loop processing solves the intermediate problem: 

Minimize: F j,] = ^ +a ( / _ J j + ^(/-jjjfaj, 

Subject to: and ^ = 0 on 5Q 



3 Returning to Figure 4, in block 127, if A> I, the result / jjfc is up scaled (2:1 ratio) by 

4 pixel replication into the new lg y the initialization for the next resolution k- 1 layer. The 

5 resolution layer is updated k = k - 1, and the algorithm proceeds by going again to block 125. 

6 If k = 1 , the result /jy, is the final output of the algorithm 1 20. 

7 By setting the parameters a = P = 0, and removing the constraint / > s, the algorithm 

8 120 reduces to the Homomorphic filtering process that was shown to be equivalent to the 

9 basic formulation of the random walk algorithms. 

10 Using the method and the system described above, the solution to Equation (1), and 

1 1 the convergence of the numerical algorithm 120 can be shown to be unique. That is, the 

1 2 variational optimization problem P, given by 

Minimize: F j /]= jj v/ f +a(/ _ g j + 
Subject to: ^ flM rf(v/ > ») = 0 on da 

1 3 a > 0 and p>0 t has a unique solution. 

14 The algorithm 1 20 has thus far been described with respect to a single channel. In 

1 5 another embodiment, the system and the method are applied to color images. When 

16 processing color images, one option is to deal with each color channel separately. Such 

17 channel-by-channel processing may be referred to as 'RGB Retinex' processing. Treating the 

18 R, G, and B channels separately usually yields a color correction effect. For example, 

19 applying RGB Retinex processing to a reddish image should modify the illumination in such 

20 a way that the red hue is removed and the resulting image is brightened and corrected. 

21 Therefore, for some image types, RGB Retinex processing improves the colors. For other 
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1 image types, such color correction may cause color artifacts that exaggerate color shifts, or 

2 reduce color saturation. 

3 In yet another embodiment, colors are mapped into a luminance/chrominance color 

4 space, such as HSV. Next, the Retinex correction is applied, but only to the luminance or V 

5 layer, and then the colors are mapped back to the RGB domain. This method may be referred 

6 to as 'HSV Retinex 5 . Color shifts in such cases are less-likely. The advantage of HSV 

7 Retinex processing is that only a single channel is processed. However, using this 

8 embodiment colors are not corrected with respect to hue (H), and saturation [s). 

9 The reflectance image R obtained by the Retinex process is sometimes an over- 

10 enhanced image. This can be explained by the facts that i) the human visual system usually 

1 1 prefers illumination in the image, and that ii) removal of all the illumination exposes noise 

12 that might exist in darker regions of the original image. 

13 In still another embodiment, a corrected version of the reconstructed illumination is 

14 added back to the constructed reflectance image. Figure 7 describes this operation. In Figure 

15 7, an algorithm 200 computes the illumination image L= exp(/) from the intensity image S = 

16 exp(.y), and the reflectance image R =S/L. This is shown where an input imageS 201 is 

17 applied to a Retinex module 202 to produce an illumination image I 203. The illumination 

1 8 image L 203 is tuned by a Gamma correction operation 204 with a free parameter y , to 

19 obtain a new illumination image U 209. A reflectance module 206 computes reflectance R 

20 207 from the input image 5201 and the illumination image L 203. The new illumination 

21 image L'209 is multiplied by the reflectance/? 207 at 210 to produce an output image S"211, 

22 where S'-L 0 R . The Gamma correction is performed by 



W 



23 where W\s the white value (equal to 255 in 8-bit images). 

24 The final result 5" is given, therefore, by: 



5' = VR = J S 

- w iLIW) S 
~~ i b -[ L IW"») 

25 For y = 1 , the whole illumination is added back, and therefore S' =S. For y = co , no 

26 illumination is returned, and 5" = R W \ which is the same reflectance image,/?, as obtained 
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1 by the algorithm 100, stretched to the interval [0, W]. This later case can also be considered 

2 as multiplying a maximum constant illumination JVto the reflectance image R. 

3 An analog to adding the illumination to the final image can be found in the 

4 homomorphic filtering approach. A linear filter for the illumination calculation in the log 

5 domain removes high-pass spatial components of s, yet also attenuates the low-pass 

6 components by a factor ofy l (where / stands for illumination). This is an analog to a gamma 

7 correction of the illumination with y = y t> since the equation for 5" can be written in the form: 

% -{wT ■* 

8 and therefore: 

s'-w = — (/ — w)+r 
Y 



= —(low - pass components) + (high - pass components). 

r * * 

9 In an ideal environment, the illumination! is assumed to be piece-wise smooth. That 

10 is, L is smooth everywhere except at strong discontinuities of the input image 5. This ideal 

1 1 can be approached in a number of ways. First, referring again to Figure (1 ), both the 

12 illumination L and the reflectance R are required to be smooth. However, these requirements 

13 contradict each other because / plus r must equal s. Thus / and r share the discontinuities. 

14 The parameter p arbitrates this contradiction and makes sure that appropriate parts of each 

15 discontinuity are apportioned to / and r. 

16 A second approach ensures that the illumination / smoothness requirement be valid 

17 everywhere except over discontinuities of s, and that the reflectance smoothness be valid over 

18 discontinuities of s. 

19 This second approach requires computing a weight function w(Vs) that obtains values 

20 close to J over most parts of the image s y and close to 0 over strong discontinuities. Equation 

21 (1) then becomes: 

22 F[/]= l^Vs^ (3) 
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1 The projected steepest descent solution to Equation (3) composed of a series of gradient steps 

2 is: 



3 



1 =/. +4pf(w(V*)p;(f l ))+l^Wv*)p5fc))-aft -*)+ 

fi'Dtfo-Y&ijMQ, -^)) + ^ Z>/((l-w(V5))D;(/ f -s 

4 and projections: 

5 li = max{li,s} 

6 Where: 

7 n is a constant 

8 D[ (0 is a forward derivative in the x-direction. 

9 D h x (.) is a backward derivative in the x-direction. 

10 D f y (.) is a forward derivative in the y-direction. 

1 1 D b y (.) is a backward derivative in the y-direction. 

12 Figure 8 illustrates an algorithm 125a, which is an alternative to the algorithm 125 in 

13 algorithm 120, that may be used to solve the functional (Equation (3)). 

14 In block 221 , the algorithm 125a calculates w(Vs), Ok and T k : 

15 "^-nfe 7 



Th 



16 where 77? is a constant threshold, 

17 and * ° _ , 

r t = r 0 -2 i+l 

18 Note that Ok changes for each resolution level-a higher ct k corresponds to a course resolution 

19 level. 

20 In block 223, the algorithm 125a calculates: 

21 G=Z>/ (w(Vs)Z) J t (l i )) + 

22 D' y (w(Vs)Z> r * (li))-a(l r s) + 

23 P D{ ((l-w(Vs)) D* (li-s)) + 

24 P I>; ((l-w(Vs)) £»; (li-s)) 
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1 In block 225, the algorithm 1 25a completes the SD iteration: 

2 I M =I 9 +M-G 

3 In block 227, the result is projected into the constraints 

4 lj = max {lj, Sk}. 

5 In block 229, the algorithm 125a checks if i = T k , and if so, processing according to 

6 the algorithm 220 ends. 
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1 In the claims: 

2 1 . An image enhancement method, comprising: 

3 capturing (121) an image; 

4 constructing (123) a multi-resolution structure comprising one or more resolution 

5 layers; 

6 processing (125) each resolution layer using an iterative algorithm having at least one 

7 iteration; 

8 projecting (149) each processed resolution layer to a subsequent resolution layer; 

9 up-calling (129) each projected resolution layer to the subsequent resolution layer; 

1 0 and using the projected resolution layers to estimate (1 30) an iflumination image. 

1 1 2. The method of claim 1, further comprising, for each of one or more iterations: 

12 calculating (143) a gradient of a penalty functional; and 

13 computing (132) an optimal line-search step size. 

14 3. The method of claim 2, wherein the penalty functional is given by: 

15 F [/]= J|V/| J + a{l-sf + /3\V(l- S f}kcfy; 

n 

16 subject to / > s and (vi, rij = 0 on dCl; wherein Q is a support of the image, cQ is an image 

17 boundary, n is a normal to the image boundary, and a and p are free non-negative real 

1 8 numbers 

19 4. The method of claim 2, wherein the penalty functional is given by: 

20 F[l] = J, (w, {Vsplf + a(l -s)+pw 2 (VjJV/ - Vsf }lxdy 

21 where wj and W2 are non-linear functions of the gradient. 

22 5. The method of claim 1, wherein the iterative algorithm is a Projected Normalized 

23 Steepest Descent algorithm. 

24 6. The method of claim 1, wherein the iterative algorithm is a Steepest Descent 

25 algorithm. 

26 7. The method of claim 1, wherein a set of constraints comprise that the illumination is 

27 greater than the image intensity, L>S. 

28 8. The method of claim 1, further comprising applying penalty terms, the penalty terms, 

29 comprising: 

30 that the illumination is spatially smooth; 

3 1 that the reflectance is maximized; 
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1 that the reflectance is piece-wise smooth. 

2 9. The method of claim 1, further comprising: 

3 computing the reflectance image based on the captured image and the estimated 

4 illumination image; 

5 computing a gamma correction factor; 

6 applying the gamma correction factor to the estimated illumination image; and 

7 multiplying the gamma-corrected illumination image and the reflectance image, 

8 thereby producing a corrected image. 

9 10. A system (1 12), embodied in a computer-readable medium, for enhancing digital 

1 0 images, 

1 1 comprising: 

12 a log module (102) that receives an input digital image S (101) and computes a 

13 logarithm s (103) of the input digital image; 

14 an illumination estimator module (104) that produces an estimate I* (107) of an 

1 5 illumination component L of the input digital image 5, wherein the estimator module (1 04) 

16 employs a construct comprising one or more resolution layers, and an iterative algorithm that 

1 7 processes each of the one or more resolution layers; and 

18 a summing node (108) that sums the logarithm s (103) and a negative of the estimate 

19 /* (107) to produce an estimate r* (109) of a logarithm of a reflectance component R of the 

20 input digital image 5, wherein a processed resolution layer is used to up-scale a subsequent 

21 resolution layer. 



18 



WO 02/089062 



PCT/US02/14987 




WO 02/089062 



PCT/US02/14987 




WO 02/089062 



3/9 



PCT/US02/14987 



DISPLAY 










IMAGE 
PROCESSOR 




IMAGE 
MEMORY 








CONVERTER 


B 









CO 

6 



PCT/US02/14987 



4/9 



INPUT 



INITIALIZE 
PYRAMID AND 
PARAMETERS 



MAIN 
LOOP 



k>l 




NO 



FINAL 
IMAGE 



FIG. 4 



WO 02/089062 



PCT/US02/14987 



5/9 



INITIALIZE 
COUNTER 
1=1 



131 



135' 



INCREMENT 
COUNTER 





FIG. 5 



WO 02/089062 



PCI7US02/14987 



6/9 



141' 



143' 



INITIALIZE 
RESOLUTION 
LAYER 



CALCULATE 
GRADIENT 
GA 



145' 



CALCULATE 
M'NSD 



147' 



COMPLETE 

NSD 
ITERATION 



149' 



125 



PROJECT 
ONTO 
CONSTRAINTS 




YES 



FIG. 6 



WO 02/089062 



7/9 



PO7US02/14987 



8^ 



GO 




6 



WO 02/089062 



PCT/US02/14987 



8/9 



221' 



223' 



INITIALIZE 
RESOLUTION 
LAYER 



CALCULATE 
GRADIENT 



225' 



COMPLETE 
SD 

ITERATION 



125a 



227 



PROJECT 
ONTO 
CONSTRAINTS 




FIG. 8 



WO 02/089062 



PCT7US02/14987 




This Page is Inserted by IFW Indexing and Scanning 
Operations and is not part of the Official Record 

BEST AVAILABLE IMAGES 

Defective images within this document are accurate representations of the original 
documents submitted by the applicant." 

Defects in the images include but are not limited to the items checked: 

□ BLACK BORDERS 

□ IMAGE CUT OFF AT TOP, BOTTOM OR SIDES 
FADED TEXT OR DRAWING 

EiBLURRED OR ILLEGIBLE TEXT OR DRAWING 

□ SKEWED/SLANTED IMAGES 

□ COLOR OR BLACK AND WHITE PHOTOGRAPHS 

□ GRAY SCALE DOCUMENTS 

□ LINES OR MARKS ON ORIGINAL DOCUMENT 

□ REFERENCE(S) OR EXHIBIT(S) SUBMITTED ARE POOR QUALITY 

□ OTHER: . 

IMAGES ARE BEST AVAILABLE COPY. 
As rescanning these documents will not correct the image 
problems checked, please do not report these problems to 
the IFW Image Problem Mailbox. 



